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Frustration has proved to give rise to an extremely rich phenomenology in both quantum and 
classical systems. The leading behavior of the system can often be described by an effective model, 
where only the lowest-energy degrees of freedom are considered. In this paper we study a system 
corresponding to the strong trimerization limit of the spin 1 /2 kagome antiferromagnet in a magnetic 
field. It has been suggested that this system can be realized experimentally by a gas of spinless 
fermions in an optical kagome lattice at 2/3 filling. We investigate the low-energy behavior of both 
the spin 1/2 quantum version and the classical limit of this system by applying various techniques. 
We study in parallel both signs of the coupling constant J since the two cases display qualitative 
differences. One of the main peculiarities of the J > case is that, at the classical level, there is 
an exponentially large manifold of lowest-energy configurations. This renders the thermodynamics 
of the system quite exotic and interesting in this case. For both cases, J > and J < 0, a finite- 
temperature phase transition with a breaking of the discrete dihedral symmetry group D% of the 
model is present. For J < 0, we find a transition temperature /\J\ — 1.566 ± 0.005, i.e., of 
order unity, as expected. We then analyze the nature of the transition in this case. While we 
find no evidence for a discontinuous transition, the interpretation as a continuous phase transition 
yields very unusual critical exponents violating the hyperscaling relation. By contrast, in the case 
J > the transition occurs at an extremely low temperature, T^~ ~ 0.0125 J. Presumably this 
low transition temperature is connected with the fact that the low-temperature ordered state of the 
system is established by an order-by-disorder mechanism in this case. 

PACS numbers: 75.10.Jm, 75.45. +j, 75.40.Mg, 75.10.Hk 



I. INTRODUCTION 

The study of frustrated quantum magnets is a fasci- 
nating subject that has stimulated many studies within 
the condensed matter community in recent yearsP^Such 
systems are assumed to be the main candidates for a rich 
variety of unconventional phases and phase transitions 
such as spin liquids and critical points with de-confined 
fractional excitations. 4 Frustration can also play an im- 
portant role in clas sical systems. The phenomenon of 
order-by-disorder^^ is the perfect example where the in- 
terplay of frustration and fluctuations produces the emer- 
gence of unexpected order. Order-by-disordcr implies 
that a certain low-temperature configuration is favored 
by its high entropy, not by its low energy. Order-by- 
disorder can also occur in a quantum system, where a 
naive argument suggests that quantum fluctuations play 
the same role as thermal fluctuations in the classical sys- 
tem, albeit there arc counterexamples where their role is 
in fact quite different.^ 

A particularly illustrative example is provided by 
the spin 1/2 antiferromagnet on the kagome lattice. 
A spin gap appe ars to be present both at zero 
magnetizatiorplSlliil (see, however, Refs. [PoT - H7|l and at 
1/3 of the saturation value where it gives rise to a 



plateau in the magnetization curve.^^T] One would be 
tempted to believe that the nature of the ground state 
is similar in both cases. However, whether the ground 
state at zero field is ordered or not is still under de- 
bate and also the existence of a plateau in the isotropic 
spin 1/2 Heisenberg model at magnetization 1/3 has been 
questioned recently! 22 ! 23 ! Nevertheless, the existence of a 
plateau at magnetization 1/3 is quite clear for easy-axis 
exchange anisotropics^^ and, using a correspondence 
with a quantum dimer model on the honeycomb lattice,^ 
the ground st ate i s identified as an ordered array of res- 
onating spins.^H 

In this paper we study an effective model that arises 
in the strong trimerization limit of the spin 1/2 kagome 
antiferromagnet! 2 ^ This model has played an important 
role in analyzin g the zero- field properties of the kagome 
antiferromagnet J2ZI25I but here we will focus on magne- 
tization 1/3 of the Heisenberg model, corresponding to 
full polarization of the physical spin degrees in the effec- 
tive model. Thus, we are left with the chirality degrees 
of freedom of the original antiferromagnet which we will 
treat as 'spin' variables. In this sense our spin system can 
be considered as a purely orbital model similar to com- 
pass models recently considered in the literature (see, 
e.g., Refs. 29 37). As an experimental realization of this 
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model a system of spin-polarized fermions trapped in a 
trimerized optical kagome lattice at 2/3 occupancy has 
been suggested!^®] j n fact, experimental realization of 
an optical kagome lattice has been reported recently,^ 
albeit using a setup which does not allow direct control 
of trimerization. 

Beyond the particular realizations of our model its 
very rich physics which results from the interplay be- 
tween classical and quantum fluctuations and frustration 
makes it an interesting model in its own right. As will 
be shown in this paper, a Hamiltonian with an (unusual) 
discrete symmetry but with a continuous degeneracy of 
the classical ground state, as it would be expected for 
a Hamiltonian with a continuous symmetry, is just one 
aspect of the rich phenomena emerging from this model. 

The present paper is organized as follows: in section 
|ll]we present the Hamiltonian and the symmetries of the 
classical and spin 1/2 cases. The Hamiltonian can be 
defined for both signs of the coupling constant J. We de- 
liberately discuss in parallel the two cases throughout the 
entire paper to point out their similarities and differences. 
The spin 1/2 case is then treated in section III by means 



(a) 



of exact diagonalization techniques, and we argue that 
a finite-temperature phase transition takes place. Since 
exact diagonalization can access only small lattices, we 
move to the classical model in section |IV| We study in 
detail the manifold of lowest-energy configurations and 
their corresponding spin-wave spectra. The effect of soft 
modes in the order-by-disorder selection mechanism is 
argued to be the origin for the phase transition of the 
J > case, in contrast to the J < case, where the 
transition is of a more conventional purely energetic ori- 
gin. In section|V]we apply Monte-Carlo techniques to the 
classical model and determine the transition temperature 
for J > and J < 0. We also analyze the universality 
class of the transition, however, only for J < since the 
transition temperature for J > turns out to be so low 
that it is difficult to access. Finally section [Vl] is devoted 
to some concluding remarks and comments. 



II. HAMILTONIAN AND SYMMETRIES 



A. Hamiltonian 



We will study the Hamiltonian 
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FIG. 1. (Color online) (a) Triangular lattice with assign- 
ment of bonds to the three different directions and underlying 
trimerized kagome lattice, (b) The two chirality states of a 
triangle, (c) Assignment of the vectors e*i to the bonds of the 
triangular lattice for the alternative representation |3| of the 
Hamiltonian. 



H = J \Y,1?Tf + £ Tfrf+^lfl* 

di,3) {{k,j}) [[k,i]] 



where 



T, = S+ + S- = 2 Sf . 



(1) 



If = u S+ + u? Sr = -Sf - V3SV , (2) 

If = J St +u>Sr = -Sf + V3SV, 

with the third root of unity w — e 21 ™/ 3 . The sums in 
|l]) run over the bonds of a triangular lattice, each cor- 



responding to one of the three distinct directions of the 
lattice, as sketched in Fig. [lja) . 

The Hamiltonian Q arises as an effective Hamiltonian 
for the trimerized kagome lattice, sketched in Fig. [ija) 
behind the triangular lattice. Our notation follows the 
derivation from the half-integer spin Heisenberg model 
for the case where the remaining magnetic degrees of free- 
dom are polarized.^ In this case, there are two pseudo- 
spin states of opposite chirality for each triangle, see 
Fig. |ljb). As reviewed in appendix [Aj plain first-order 
perturbation theory of the S z -S z interactions between 



3 



two triangles yields Eq. ([I]) where the exchange constant 
J is proportional to the inter-triangle exchange constant 
of the kagome lattice and would thus typically assumed 
to be antiferromagnetic (J > 0). Note that we have cho- 
sen a convenient normalization of J. A similar derivation 
starting from a Fermi gas with two atoms per trimer also 
leads to the Hamiltonian 0P^ 

Due to the two possible chiralities on each triangle, the 
pseudo-spin operators Si should be co nsider ed as quan- 
tum spin- 1/2 operators. The derivation j 26 l 38 l also suggest 
a positive J > to be more natural. In this paper we will 
relax these constraints and, for reasons that will become 
clear later, consider also classical spins, i.e., unit vectors 
Si, and the case J < 0. 

Note that the Hamiltonian is not symmetric un- 
der reflections of the lattice. Our conventions agree with 
those of Ref. [26] where this Ha milton ian appeared first, 
while some more recent workJ^ H 40 * 42 ! use a reflected con- 
vention for the chirality. Note furthermore that our con- 
ventions for J differ by a factor 4 from previous studies 
of the model 0P^ 

It may also be useful to represent the Hamiltonian 
in a more compact form 4 ^ 

H = 4jJ2(s i -§ i ) , (3) 

(id) 

where the vectors e*j are indicated in Fig. [ljc) : for each 
bond one has to choose e*j and ej as in the corresponding 
bond of the bold triangle. For example, for each horizon- 
tal bond (i, j) , one needs to choose e*j = ( J j for the left 

site i and &j = \ ^ ^J^J ^ or ^ e rl S^ site j- 

Models which are very similar to have recently been 
studied in the context of spin-orbital models (see, e.g., 
Refs. EM37]). 

B. Symmetries 

The Hamiltonian has the following symmetries on 
an infinite lattice: 

1. Translations T x , T y along the two fundamental di- 
rections of the lattice. 

2. Simultaneous rotation i?27r/3 of the lattice and all 
spins around the z-axis by angles of 27r/3 (the latter 
rotation amounts to a cyclic exchange of Tf-, Tf 3 , 
and if). 

3. A rotation by 7r around the z-axis in spin space: 
P : Sf i-> -Sf, Sf i— > -S? while keeping the 
lattice fixed. 

4. A spatial reflection combined with rotation of all 
spins around a suitable axis in the x-y-plane by 
an angle tt. One particular choice is I : Sf 



Sf, Sf h> -Sf, Sf h> -Sf, combined with a 
reflection of the lattice along the dashed line in 
Fig- [3a). 

5. For spin 1/2, there is another symmetry imple- 
mented by 

Q=ri( 2 ^)- w 

i 

Conservation of Q means that the number of spins 
pointing up (or down) along the z-axis is a good 
quantum number modulo two. This conservation 
law is most easily verified by observing that the 
interaction terms in always invert a pair of spins 
in an eigenbasis of S z . 

The choice of factors in ensures that Q 2 = 1. Fur- 
thermore, one has P% w / 3 = P 2 = I 2 = 1. i?27r/3 and P 
together generate the abelian group Zg = Z3 x Z2, as 
described for instance in Ref. [351 The combination of 
^2u-/3i P, and I generates the dihedral group Dq, which 
is non-abelian (IR 2lT /3l = ^2^/3)- Finally, -Raw/a and 
/ generate the symmetric group £3 which can be traced 
to the point-group symmetry of the underlying kagome 
lattice. The operators i?2u-/3j P, and I leave the Hamil- 
tonian Q invariant irrespective of the value of the spin 
quantum number. Thus, the group D e is a symmetry 
also of the classical variant of the model. 

The symmetries P and Q are not present in the un- 
derlying kagome lattice, hence they should be specific to 
the lowest-order approximationPSEH Indeed, at least in 
the derivation from the Hciscnbcrg model one observes 
that already the next correction^ breaks the symmetries 
P and Q. 

Now let us consider the consequences of the combina- 
tion of I and Q for the spin- 1/2 model on a finite lattice 
with N sites. Then the relation I S? I = —S? leads to 

IQI = (-l) N Q. (5) 

Since the eigenvalues of Q are q = ±1, this implies that 
/ is an isomorphism between the subspace with q = — 1 
and the subspace with q = 1 for odd N and spin 1/2. 

III. QUANTUM SYSTEM: EXACT 
DIAGONALIZATION FOR SPIN 1/2 

In this section we will present numerical results for 
the Hamiltonian with spin 1/2. We impose periodic 
boundary conditions and use the translational symme- 
tries T x , T y in order to classify the states by a momentum 
quantum number k. We only consider lattices which do 
not frustrate a potential three-sublattice order, i.e., only 
values of N that are multiples of three. For the system 
sizes considered already in Ref s . l39l and 140} we will use 
the same lattices. In particular, the N — 12, 21, and 24 
lattices are shown in Fig. 9 of Ref. HOI Furthermore, we 
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will consider the N = 27 lattice which can be found, e.g., 
in Fig. 1 of Ref. HD 

Let us briefly discuss the consequences of the other 
symmetries mentioned above. We did not make explicit 
use of P, although it is present for any lattice. However, 
the symmetry Q (which is also present for any lattice, see 
Q) is easily implemented if we work in an S^-eigenbasis. 
Concerning the rotation i?27r/3i it is not possible to find 
lattices for all N such that it is a symmetry. If i?27r/3 
is a symmetry, we use it to select one representative k 
for all equivalent momenta. Finally, the presence of the 
symmetry / is more delicate. We have performed com- 
puter checks and found that most of the lattices under 
consideration have a suitable spatial reflection symme- 
try, ensuring that / is a symmetry The only exception 
is the N = 21 lattice where there is no such reflection 
symmetry. Nevertheless, we find the same spectra in the 
subspaces with q = — 1 and q = 1 also for N = 21. 
Therefore, for N odd we can choose representatives for 
all symmetry sectors in the subspace with q = 1. 

For N < 21 the translational symmetries and Q lead 
to matrices with dimension up to 49 940 and we can ob- 
tain all eigenvalues. Dimensions increase up to 2 485 592 
for N — 27. In this case, we have used the Lanczos 
method to compute the n lowest eigenvalues in each sec- 
tor. Mainly for reasons of CPU time, we restrict to 
n re 70 (150) for N = 27 (24) and J > 0. 



A. Low-lying spectra 

Let us first look at the spectra. In order to take de- 
generacies into account, in Fig. [2]we show the integrated 
density of states, i.e., the number of states with energy 
less or equal than AE above the ground state. 

Panel (a) of Fig. [2] shows results for J < 0. These 
results extend previously presented results^ for N = 12, 
18, and 21 to higher energies and larger N. One observes 
that there are at most 8 states for energies AE < 3 | J\ 
with a substantial density of states setting in at higher 
energies. This suggests a thermodynamic gap ps 3 \J\. 

Fig. [2|b) shows the density of states for J > 0, ex- 
te nding previously published results for N — 18, 21, and 
24j 39 l 4 ° l In this case, we observe a large density of states 
at substantially lower energies than for J < 0. This large 
density of states is reminiscent of the large density of non- 
magnetic excitations observed in the spin-1/2 Heisenberg 
antiferromagnet on the kagome lattice, both at zero mag- 
netic fielcP™ and on the one-third plateau.^ In partic- 
ular the N = 27 data presented in Fig. [2jb) shows a 
large density of states for AE > 0.02 J. On the other 
hand, one observes at most 8 states with AE < 0.012 J 
in Fig. ^b) for a given system size N. From these ob- 
servations we infer that a gap is at most on the order of 
0.02 J if present at all. 

Since an ordered ground state breaks the symmetry 
group Dq, such a ground state should be six- fold degen- 
erate. Indeed, classical and semiclassical considerations 
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FIG. 2. (Color online) Integrated density of states, i.e., num- 
ber of states with energy less or equal than AE above the 
ground state, for (a) J < and (b) J > 0. 



predict a six-fold degeneracy in an ordered state (see sec- 
tion 



IV below). However, there is no clear separation of 



6 low-lying states from the remainder of the spectrum 
for J < (see Fig. [2^a)), and even less so for J > 
(Fig.[2][b)). The considered lattice sizes may be too small 
to observe the expected low-energy structure of the spec- 
trum. However, correlation functions exhibit pro noun ced 
120° correlations already on these small lattices J22HS2I 



B. Specific heat 

The specific heat C can be expressed in terms of the 
eigenvalues of the Hamiltonian. Since we have all eigen- 
values for N < 21, it is straightforward to obtain the 
specific heat for all temperatures and both signs of J. 
Fig.[3]shows the results of the specific heat per site C/N. 
The case J < is shown in panel (a). There is a finite- 
size maximum slightly above T ss |J|. The large finite- 
size effects which are still observed here are consistent 
with a phase transition around T « |J| in which case 
C should become non- analytic for N — > oo. Because of 
a possible phase transition, we have tried to obtain a 
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low-temperature approximation to the specific heat for 
N > 21 , J < by keeping low-energy states. However, 
for N = 24 even 12 462 low-lying states going up to ener- 
gies as high as AE < 12.6 \J\ turned out to yield a specific 
heat which has sufficiently small truncation errors only 
for temperatures T < 0.9 |J|. This result for N = 24 
(also included in Fig. pjta)) clearly does not include the 
maximum of the specific heat C. 

Now we turn to the case J > which is shown in panel 
(b) of Fig. [3j In this case, finite-size convergence at high 
temperatures is much faster than for J < 0. This fast 
convergence indicates that there is no phase transition 
associated to the high-temperature maximum (the posi- 
tion of this maximum is at T rj 2.1225 J with a value 
C « 0.105717iV for N = 21). The reduced finite-size 
effects and the smaller value of C at the maximum re- 
flect that there is a substantially smaller energy scale for 
J > as compared to J < 0. 

For J > 0, a second peak emerges in the specific heat at 
low temperatures, see Fig. [3^b). In order to investigate 
this in more detail, we use again the low-temperature 
approximation for the specific heat obtained from the 
low-lying part of the spectrum. For N = 24 and 27, we 



have used a total of 7 029 and 3 906 eigenvalues, respec- 
tively (the N — 24 data is included in Fig. [3jb), but it 
is difficult to see there since it is valid only at very low 
temperatures). Fig. pi shows the specific heat divided by 
temperature (panel (a)) and the entropy per site (panel 
(b)) in the low-temperature region for J > and system 
sizes N — 18, 21, 24, and 27. Our result for the spe- 
cific heat with N = 21 obtained from the full spectrum 
agrees with a previous result for iV = 21 based on ap- 
proximately 2 000 low-lying statesPl The finite T = 
limit of the entropy for N = 21 and 27 in Fig.[4|b) corre- 
sponds to the two-fold degeneracy of the ground state for 
these system sizes, see Fig. 2jh). Although the maximum 
value oiC/T increases with increasing N, there are non- 
systematic finite-size corrections to the position of this 
maximum. Thus, we can only conclude that if there is a 
finite-temperature ordering transition for J > 0, it should 
have a very low transition temperature T c < J/100. 

Fig. Kb) shows that there is a remarkably large en- 
tropy S/N w 0.2 ... 0.3 associated to the finite-size low- 
temperature peak of the specific heat. This is comparable 
to the entropy associated to the degeneracy of the classi- 



cal ground states, see section IV C below. Therefore, this 
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observation lends further support to the interpretation^! 
of the low-energy states for S = 1/2 in terms of the clas- 
sical ground states for J > 0. 



IV. CLASSICAL SYSTEM: LOWEST-ENERGY 
CONFIGURATIONS AND SPIN-WAVE 
ANALYSIS 

We will now proceed with a discussion of the low- 
energy, low-temperature properties of the classical vari- 
ant of the model |l]), i.e., we will assume that the Si are 
unit vectors. We will parametrize the spin at site i by 
angles oti and 7$: 



cos 7, cos a, 
cos 7i sin a.i 
shi7j 



(6) 



Because the z-components do not contribute to the en- 
ergy, configurations with extremal energy should have 
spins lying in the x-y-plane (7$ = 0). The energy E({a 
for a given set of angles {cti}, 7» = is obtained from (|1 
by identifying 



= 2 cos (a,) , 
Xf = 2 cos (a, + fi) 
Tf = 2 cos (an - fi) 



(7) 



with fi = 2tt/3. 

We will further be interested in small fluctuations 
{a.i + £.;}, {ji — ?j} around a ground-state configuration 
{«i,7i = 0}. The energy can then be expanded as 



E ({a, + q}) = E ({a,}) + ^ a M hJ e, 
+ E zz + 0({e i ,£ j } 3 ) . 



(8) 



Here, E zz is a diagonal quadratic function of the out-of- 
plane fluctuations it which, to quadratic order, decouples 
from the relevant degrees of freedom e,. The eigenvalues 
fi of the symmetric matrix Mij correspond to the spin- 
wave modes. The fact that {o^} describes a ground state 
implies /j > 0. We will call a mode with /j = 'pseudo- 
Goldstone mode'. 



A. Ground states with small unit cell for J < 

Let us first consider the case J < 0. Then a ground 
state is given by a certain three-sublattice configura- 
tion where the angles b etwee n different sublattices 
differ by multiples of 27r/3pMD pjg j^J g^ows sucn a 
low-temperature configuration as a snapshot which was 
taken during a Monte-Carlo simulation (details to be 
given in section [V] below). The energy of such states 
E^ = 6 J N is invariant under global rotations of the 




FIG. 5. (Color online) Snapshot of a configuration during a 
simulation for J < at T = 1CT 3 j Jj on a 12 x 12 lattice. 
Periodic boundary conditions are imposed at the edges. 



spin configuration in the x-y-plane, i.e., there is a one- 
parameter family of ground states (note that this invari- 
ance under a continuous group is not a symmetry of the 
Hamiltonian). We parametrize this global rotational de- 
gree of freedom by an angle a of the spins on one sublat- 
tice. Using a three-site unit cell, we can exploit invariancc 
of this ground state under translations to represent the 
matrix pi in Fourier space by the following 3x3 matrix 




where 



( 6 

M = J 

\-2B* 



A = c ikl sin 2 (a) + e l k2 sin 2 (a + fi) 
+e' k3 sin 2 (a -CI) , 



(9) 



B = 



-i k± 



sin z (a - n)+c- lk2 sin^ (a) 



and 

fci = 



+c- ik3 sin 2 (a + fi) , 
C = c ikl sin 2 (a + fi) + e l fea sin 2 (a - fi) 
+e 4fc3 sin 2 (a) , 



(10) 



Let us analyze now the effect of the fluctuations by 
computing the free energy associated with (|8j). To this 
end we can compute the partition function 



-PH z 



\{de{k) 



(12) 
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FIG. 6. (Color online) Low-temperature contribution to the 
free energy ( 15 1 for J < of the fluctuations above the 120° 



ground state as a function of the spin angle a. 



where f"(k) are the eigenvalues of (|9j) and Z zz is the 
Gaussian integral over the N quadratic variables corre- 
sponding to the out-of-plane fluctuations. 
Performing the Gaussian integral we get 



-ma 



which yields the free energy as 
Ann/3 1 



F = H Q + F Z 



2(3 2(3 



i , k 



(13) 



(14) 



where F zz = ~\nZ zz /(3. 

The low-temperature behavior is therefore determined 
by the following contribution of the fluctuations to the 
free energy 

= £ In = £ In detM, (15) 



where M is the matrix ( 9 1 
gral ( 15 ) is shown in Fig. 6 



The result of the inte- 
J 7< (a) is a 2 7r/3-periodic 
function since the spin angles on the different sublat- 
tices differ by 2tt/3. Hence, it is sufficient to consider 
a e [0, 2 it/ 3}. One observes that J r< and thus the low- 
temperature limit ( 14 1 of the free energy has minima at 
a = (2 n + 1) 7r/6, n = 0, 1, . . . , 5. This implies that the 
120° classical ground-state configuration locks in at these 
angles for T — > 0. This lock-in can indeed be verified 
in histograms of Monte-Carlo simulations, see Ref. 03] 
for a planar variant of this model and also Fig. [TT] be- 
low. Lock-in of the classical ground-state configuration 
at a = (2n+l)7r/6 follows also from the semiclassical 
approach.^ In this approach the ground-state energy is 
given by 



^semclass.H - E , 



< 

class. 



6JS + YN u f@,a).. 



\j\ 




2 7r/(3 V3) 



FIG. 7. (Color online) The three eigenmodes /» for the 120° 
ground state with J < 0, assuming a lock-in of the ground 
state at a — (2n + l)n/6. Note that the shaded surfaces 
extend only over the first Brillouin zone. 



where ujf (k, a) are the three sheets of spin-wave (SW) 
frequencies obtained from the linear Holstein-Primakoff 
approximation and where the sum over k runs over 
the magnetic Brillouin zone. The SW frequencies 
are connected with the classical eigenmodes f^ ass t by 



icLs.,» = K < (fc,a)) 2 /(245 2 |J|). One finds that 

-^stmciass. ( a ) nas mm i m a at a = (2n + 1)tt/6. The 
expressions for /<j ass 8 -(a) are too cumbersome to be 
given explicitly, instead Fig. [7] shows a plot of the three 
cigenfrequencies fi at a = (2n + l)ir/6. The lowest 
sheet has a unique quadratic minimum at k = with 
/cis,.i(fc)«9|^|fc 2 /8foi'fc«0. 



(16) 



B. Ground states with small unit cell for J > 

N ow we turn to the case J > 0. There is a first ground 
stat d 39 * 4 ^ described by on — a with an arbitrary angle a. 
This 'ferromagnetic' state has energy E^H° — — 3 J N. 
However, for J > there is another ground state with a 
small unit cell j 39 l 4 ° l again with three sublattices where the 
angles on between different sublattices differ by multiples 
of 2 7r/3. Also the energy = —3 J N is invariant 

under global rotations. The latter state is illustrated by 
the global structure of Fig. [8] which shows a snapshot of 
a low-temperature configuration taken during a Monte- 
Carlo simulation (details again to be given in section |v| 
below). Note that the sense of orientation around a tri- 
angle, i.e., the chirality of the spins in Fig. [8] is exactly 




FIG. 8. (Color online) Snapshot of a low-temperature con- 



figuration during a simulation for J > at T 



10" 



a 12 x 12 lattice. Periodic boundary conditions are imposed 
at the edges. Different colors are used for each of the three 
sublattices. 



the opposite of Fig. [5] 

The ferromagnetic state is the simplest case for the 
computation of fluctuations since Mij is diagonalized by 
a Fourier transformation. One finds the modes 



f ferro I jL jL \ 
3 class. V ft tt> h y) 

4 J 



= - + sin (a) sin (a — 12) cos (fci) 

+ sin (a + 12) sin (a — 12) cos (fe) 
+ sin (a) sin (a + 12) cos (^3) , (17) 

with the ki defined in Eq. ( [TT| ). As for the case J < 
0, the classical frequencies f^^. are proportional to 
the squares of the SW frequencies u/ erro obtained from 
a linear Holstein-Primakoff approximation®! f*fj£° = 

(w fcrro ) 2 /(l2^ 2 J). 

By computing the contribution of the modes f^f™ 
to the free energy we find minima for a = nir/3, n = 
0, 1, 2, . . ., so that the spins in the ferromagnetic struc- 
ture lock in to the lattice directions of the triangular lat- 
tice. For the lock-in values of a, f^H° depends only on 
one of the ki and has a line of zeros in the perpendicular 
direction in momentum space. 

The three-sublattice state leads to the following 3x3 
matrix in Fourier space: 



+ (e lkl sin(a + 12) + e ik3 sin(a - 12)) sin(a) , 



( e sm 
B = c- lkl sm(a + 12) sin(a - 12) 



c 

+ (e~ zk3 sin(a + 12) + e' tk2 sin(a - 12)) sin(a) , 
C = c lk3 sin(a + 12) sin(a - 12) (19) 
+ (c ife2 sin(a + 12) + e lkl sin(a - 12)) sin(a) . 



For a — nir/3, diagonalization of (18) leads to three 
completely flat branches 



1 (k) — , /class 2 



(k) 



class., 3 



{k) = lj. (20) 



•I cla 

In particular the lowest branch f^ ass 1 = corresponds 
to a branch of soft modes. In real space these soft modes 
correspond to the rigid rotation of one single triangle.^ 
Note that there is no such flat branch of soft modes for 
a value of a which is not an integer multiple of 7r/3. 

When computing the contribution of fluctuations 
around these configurations (a = nir/3) to the free en- 
ergy, one finds that one third of the modes are quartic 
instead of quadratic. This yields a free energy of the 
form: 



F = H Q + F Z 



F T 



(21) 



where, again, F zz ~ N In (3/(2 f3) corresponds to the triv- 
ial contribution of out-of-plane fluctuations and (compare 
also Ref. H]) 



Nlnp N In (3 



3/3 



12/3 



+ 



(22) 



At low temperatures, this term dominates the free energy. 
The flat branch of soft modes reduces the coefficient of 
In /3//3 from N/ 2 as in the case of only quadratic modes 
(compare §L4$) to AT/3 + N/12 = 5 N/12. This implies 
two things: firstly, the angles of the 120° state should 
lock in at a = nir/3 for low temperatures. Secondly, 
a thermal order-by-disorder mechanism should favor the 
120° state over the ferromagnetic state for T — >• 0. 

As in the case of the ferromagnetic state one finds 
that the relation f^ laas i = (w>) 2 / (12 S 2 J), where w>, 
i = 1, 2, 3, are the SW frequencies o btaine d from a lin- 
ear Holstein-Primakoff approximationpSESI holds for ar- 
bitrary values of a. Using the results for u> °(k, a) and 
u}f"(k,a) to calculate semiclassical ground-state energies 
^SdfsJ") and E> mcUss (a) in the same manner as 
in ( |16[ ) one finds that both are minimal at a = mr/3 
and that £> mclas Jn7r/3) < i& ass >7r/3). Thus the 
semiclassical approach is fully consistent with the classi- 
cal findings. 



3 2A 2B 
V/ = J I 2 A* 3 2C 
s 2B* 2C* 3 

where 

A = e l k2 sin(a + 12) sin(a - 12) 



(18) 

C. Enumeration of ground states for small N 

Direct computer inspection of all states with a ngles 
a,, = n,7r/3, n t € {0,1,2,3,4,5} for N = 12 showsPMI 
that there are further ground states for J > 0. On the one 



9 



TABLE I. Number T>n of classical ground states for J > 
on a lattice of size N with one angle fixed. The n g in the 



decomposition T>n 



denote the number of ground 



states with g pseudo-Goldstone modes. 



N V h 



12 40 = 6i + 31 2 + 2 3 + U 0.3074 

15 102 = 60i + 20 2 + 20 3 + 2 5 0.3083 

18 286 = 92i + 112 2 + 51 3 + 30 4 + U 0.3142 

21 688 = 260i + 2IO2 + 203 3 + 14 5 + 1 7 0.3111 

24 1838 = 384i + 958 2 + 199 3 + 280 4 + 16 6 + Is 0.3132 

27 5054 = IO681 + 972 2 + 2257 3 + 351 4 + 378 5 0.3158 
+ 9 6 + I87 + lg 



hand, CPU time for a similar enumeration of all 6^ such 
configurations becomes prohibitively big for N £ 15. On 
the other hand, all known ground states turn out to have 
mutual angles which are multiples of 2 tt/3. Furthermore, 
we eliminate a global rotational degree of freedom by 
fixing one arbitrary angle «o = tt- Then there remain 
just S^ -1 configurations with a* £ {7r/3,7r, —tt/3} to be 
enumerated. Direct enumeration of these 3 JV ~ 1 config- 
urations can be carried out with reasonable CPU time 
for N < 27, but becomes quickly impossible for larger 
N. We have therefore performed such enumerations for 
N < 27, using the same lattices as in section |TlT| 

The number of ground states XV determined in this 
manner is given m Table [I] for J > 0. Note tha t the 
ordered states which we have described in section IIVBI 
are just two of the 2?jv states, but there are many fur- 
ther ground states which can be interpreted as de fects in 
and domain walls between the ordered statesP^^Indeed, 
also closer inspection of the snapshot shown in Fig. [8] re- 
veals the presence of defects in the three-sublattice struc- 
ture. The 240 = 6T>i 2 states described previousl} * 39 * 40 ! 
for N = 12 are recovered by a global rotation of the an- 
gles such that ctQ takes on the six values ao = mr/3. 
The last column of Table [[] lists In (V N ) /N. The fact 
that these numbers stay almost constant indicates a fi- 
nite ground-state entropy per site slightly above 0.3 in 
the thermodynamic limit. 

It is straightforward to derive the N x N matrix Mij 
defined in Q for any ground state and diagonalize it. 
Among the eigenmodes fi, one can then identify the g 
pseudo-Goldstone modes fi — and in turn count the 
number n g of ground states with g pseudo-Goldstone 
modes. These numbers are also given in Table [I] in the 
form T>n = 'YligTig. One observes that all ground states 
have at least one pseudo-Goldstone mode. There are at 
most N/3 pseudo-Goldstone modes, and there is only 
one ground state with this maximal number of pseudo- 
Goldstone modes corresponding to the three-sublattice 
state described in section IV B (apart from N — 15; how- 
ever, this lattice is special in that it has a period three 
translational symmetry T 3 = 1). 

There are iV/3 ground states which differ from the per- 



fect three-sublattice state by a rigid rotation of the spins 
on certain triangles by an angle 2 tt/3, and another A^/3 
ground states where the spins on a different set of trian- 
gles are rotated by — 2 7r/3. These ground states have two 
pseudo-Goldstone modes less than the three-sublattice 
state. The data in Table [I] show that these 2N/3 config- 
urations with N/3 — 2 pseudo-Goldstone modes account 
for all states with the second largest number of pseudo- 
Goldstone modes for N > 21. This indicates that ground 
states deviating from the homogeneous three-sublattice 
state are obtained at the expense of pseudo-Goldstone 
modes. At finite temperature such "inhomogeneous" 
ground states are then penalized by an entropic cost be- 
cause of the reduced number of pseudo-Goldstone modes. 
This indicates that ground states with bigger deviations 
from the three-sublattice ground state have a higher free 
energy for small T since they have less soft modes. Thus, 
the global minimum of the free energy is the perfectly 
ordered state with many close-by configurations which 
deviate only locally from the perfect order. These ar- 
guments predict a thermal order-by-disorder selection of 
the 120° state among the macroscopic number of ground 
states for T -> and J > 0. 

The above enumeration procedure can also be per- 
formed for J < 0. In fact, in this case we have carried 
it out twice, first with on £ {tt/3, it, — 7r/3} and o?o = t, 
and then again with all a<j shifted by it/ 6 in order to 
match the lock-in predicted in section |IV A| In sharp 
contrast to the large degeneracy found for J > 0, we 
confirm that the ground state is unique (up to global ro- 
tations of all angles) for J < and thus identical to the 
three-sublattice state described in section |IV A| Diago- 
nalization of the corresponding N x N matrix M yields 
one pseudo-Goldstone mode, in complete agreement with 
the fi shown in Fig. [7] which have just one zero, namely 

icL,,i(£ = o) = o- 



V. CLASSICAL SYSTEM: MONTE CARLO 
ANALYSIS 

Section |IV| already contained some discussion of the 
low-temperature properties in the classical limit. The re- 



sults of section III lead us to suspect a finite-temperature 
transition in the quantum system at least for J < 0. In- 
deed, a finite-temperature phase transition is allowecP^ 
since the model has only discrete and no continuous 
symmetries (compare section |n]). Since such a finite- 
temperature phase transition should be a classical phase 
transition, we may hope to gain insight into its uni- 
versality class by studying the classical counterpart of 
the model. This motivates us to present results of a 
Monte-Carlo simulation of the classical model, treating 
the spins Si as classical 0(3) vectors. Simulations were 
performed on square lattices with diagonal bonds (topo- 
logically equivalent to the triangular lattice) and periodic 
boundary conditions. 

First, we have used a standard single-spin flip 
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Metropolis algorithm (see, e.g., Ref. |4"B"|) . Some results 
obtained from such simulations have already been pub- 
lished in Ref. 031 but the results to be presented here 
have been obtained from new runs using the 'Mersenne 
Twister' random number generatorPH In order to deter- 
mine error bars, we have used between 100 and 400 in- 
dependent simulations for J < 0. 

For J > 0, the standard single-spin flip algorithm turns 
out to be no longer ergodic for temperatures T < 0.1 J. 
Such problems are in fact expected in view of the large 
ground-state degeneracy which we discussed in section 



IV C In this region we have therefore used the parallel 
tempering Monte Carlo method (also known as exchange 
Monte Carlo - see Refs. |4"8H5T1 and references therein). 
In this framework, n simulations are performed in paral- 
lel, each at a different temperature using the standard 
Metropolis algorithm. Periodically, the exchange be- 
tween the configurations of two simulations consecutive 
in temperature is proposed and accepted depending on 
the energy balance of such a move. A careful choice of 
the temperature points allows each configuration to shuf- 
fle through the entire temperature range during the sim- 
ulation, greatly reducing the probability of getting stuck 
in a local minimum of the free energy. In principle, this 
allows an efficient exploration of the phase space, while 
not having to wait for rethermalization of the systems af- 
ter each configuration switch. Strategies to optimize the 
choice of the temperature grid have been proposed (see, 
e.g., Ref. I5ip . but we simply opted for constructing a 
fine-grained temperature set using the rule of thumb that 
the probability for a configuration switch to be accepted 
should always be at least around 70% to 80%. The result- 
ing grid consists of at least 96 points for T/j£ [0.01, 0.7], 
where of course most points lie in the low-temperature 
region. After an initial thermalization, observables are 
sampled until convergence of their error bars is observed 
(the typical duration of a simulation being at least 3 x 10 7 
Monte-Carlo sweeps per system). We would like to insist 
that even if parallel tempering is adequate to the task of 
studying the low-temperature properties of such a highly 
degenerate frustrated magnet, it is still by no means easy 
to obtain physically relevant data at such a low temper- 
ature for continuous spherical spins as we shall see later. 



A. Specific heat 

Fig. [9] shows results for the specific heat of the classi- 
cal model for J < (panel (a)) and J > (panel (b)). 
Statistical errors should be at most on the order of the 
width of the lines. Although all lattice sizes are bigger 
than those used previously for the quantum model, there 
are remarkable similarities of the specific heat of the clas- 
sical model shown in Fig. [9] with the specific heat of the 
quantum model, see Fig. [3} For J < 0, a singularity 
seems to develop in the specific heat for temperatures 
around T sa 1.5 |J|, signaling a phase transition. For 
J > 0, there is also a broad maximum at 'high' temper- 
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FIG. 9. (Color online) Specific heat per site for the classical 
model with (a) J < and (b) J > 0. Error bars do not exceed 
the width of the lines. 



atures T ~ 0.3 J. The finite-size effects for the latter 
maximum are small indicating that this does not corre- 
spond to a phase transition. In this case, the interesting 
features of the specific heat lie in the low-temperature 
region, as displayed in Fig. 10 ^b). As the system size 
increases, one can see that a small peak builds up in the 
specific heat for T « 0.02 J. This seems to indicate that 
a phase transition might occur around that temperature, 
two orders of magnitudes smaller than for J < 0. We 
are unfortunately not on a par with the J < data, as 
the CPU requirements are too steep to secure relevant 
data for systems larger than 27 x 27 sites even though 
the specific heat is a comparably robust quantity, and it 
is clear that other observables are needed to conclude on 
the existence of this phase transition. 

An important difference between the S = 1/2 and the 
classical model arises at low temperatures: the specific 
heat of the quantum system has to vanish upon approach- 
ing the zero temperature linvr->o C/N — 0, while due to 
the remaining continuous degrees of freedom, the specific 
heat of the classical system approaches a finite value for 
T — > 0. For J < 0, the equipartition theorem predicts 
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FIG. 10. (Color online) Zooms of the specific heat per site for 
the classical model: (a) in the vicinity of the maximum for 
J < 0, (b) in the low-temperature region for J > 0. Statistical 
errors do not exceed the width of the lines or the size of the 
symbols, respectively. 



N=6x6 

N=12x12 

A/=18x18 




FIG. 11. (Color online) Histogram of in-plane angles <j>i for 
J < 0. Averaging has been performed over 1000 independent 
configurations at T — 10~ 3 \J\. 



90 x 90 curves in Fig. [lO^a) demonstrate that the maxi- 
mum value of the specific heat starts to decrease as one 
goes to system sizes beyond N = 45 x 45. This implies 
that the exponent a which characterizes the divergence 
of the specific heat at the critical temperature is very 
small or maybe even negative. 



B. Sublattice order parameter, Binder cumulant, 
and transition temperature for J < 



an N/2 contribution to the specific heat for each trans- 
verse degree of freedom which yields limr_j.o C/N = 1, 
in excellent agreement with the results depicted in the 
panel (a) of Fig. [9j For J > 0, one must take into ac- 
count the fact that the flat soft-mode branch of the three- 
sublattice state is expected to contribute only N/12 to 
the specific heat. Thus for J > one should expect 
limiM-o CyiV = H/12 = 0.916666.... As can be seen in 
Fig. [Toj^b) , we observe a specific heat lower than one in 
the low-temperature region along with a downward trend 
as T goes to for all the system sizes studied. However, 
according to the data which we have at our disposal, it 
seems that one would have to go to very low tempera- 
tures T < 10 -3 J in order to verify the prediction for the 
zero-temperature limit. 

Returning to the finite-temperature transition for J < 
0, Fig. [lO^a) shows a zoom into the relevant tempera- 
ture range, including data for up to N = 90 x 90 spins. 
At these bigger system sizes, the position of the maxi- 
mum continues to shift to lower temperatures and the 
maximum sharpens. However, the N — 45 x 45 and 



According to subsection |IV A| we expect that the phase 
transition observed for J < is a transition into a three- 
sublattice ordered state. This ordering is indeed exhib- 
ited at least at a qualitative level by snapshots of Monte- 
Carlo simulations at low temperatures (compare Fig. [5| . 
In addition, one observes in Fig.[5]that the spins are lying 
essentially in the x — y-plane for low temperatures. 

Furthermore, we expect a lock-in of the spins to one 
of 6 symmetrically distributed directions in the plane at 
low temperatures (compare Fig. [6| . The latter prediction 
is indeed verified by the histogram of the angles of the 
in-plane component of the spins fa at low temperatures 
shown in Fig. [IT] Note that the histogram is rather flat 
for the smaller lattices (in particular the N — 6 x 6 lat- 
tice) and sharpens noticeably as the lattice size increase 
to N = 36 x 36 (the largest lattice which we have con- 
sidered in this context). The fact that the lock-in occurs 
only on large lattices can be attributed to the replace- 
ment of the sum over k in (151 by an integral being a 



good approximation only for large lattices. 

To test for the expected three-sublattice order, we in- 
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FIG. 12. (Color online) Square of the sublattice order param- 
eter Trig = ^A? s 2 ^ f° r the classical model with J < 0. Error 
bars do not exceed the width of the lines. 

troduce the sublattice order parameter 



(23) 



iec 



where the sum runs over one of the three sublattices C 
of the triangular lattice. Fig. [12] shows the behavior of 
the square of this sublattice order parameter for J < 0. 
One observes that the sublattice order parameter indeed 
increases for T < 2 |J| and goes indeed to m 2 = 1 for 
T — s- 0, as expected for a three-sublattice ordered state. 
Inclusion of larger lattices (up to N = 90 x 90) allows one 
to restrict the ordered phase to T < 1.7 \J\. However, 
more accurate estimates for the transition temperature 
can be obtained in a different manner. 

A useful quantity to determine the transition into an 
ordered state accurately is the 'Binder' cumulant! 46 l 52 l 53 l 
associated to the order parameter ( 23 1 via 



U s = 1 



3 (M s 4 
5 / M. ? 



(24) 



We have chosen the prefactor in (24 1 such that U s = for 
a Gaussian distribution around zero of the order param- 
eter P(M S ) = (f) 3/2 exp (-cNlfj. Such a distribution 

is expected at high temperatures, and we expect U s — > 
for T 3> | J\. Conversely, in a perfectly ordered state one 

will have ^M s 4 ^} = ^M s 2 ^> such that U s = 2/5 in this 

case. Hence, for an ordered state we expect U s ps 0.4 for 
T<T C . 

Fig. [13] shows results from classical Monte-Carlo sim- 
ulations for the Binder cumulant U s of the system with 
J < 0. First, the broad temperature range shown in 
Fig. [l3^a) confirms that indeed U s w 0.4 in the ordered 
low-temperature phase and U s ~ for high temperatures. 
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FIG. 13. (Color online) Binder cumulant U s for the classical 
model with J < 0: (a) global behavior, (b) in the vicinity of 
the critical temperature and for bigger lattices. Error bars in 
panel (a) do not exceed the width of the lines. 



The transition temperature can now be accurately ex- 
tracted from t he cross ings of the Binder cum ula nts at dif- 
ferent sizes j\r ! 46 l 52 l 5 3 For this purpose, Fig. 13 ^b) zooms 
in to the relevant temperature range, including bigger 



system sizes N. Although the crossings between any pair 
of system sizes Ni and N 2 fall into a narrow temperature 
window, there still remains a small residual dependence 
on the sizes N\ and N 2 considered. In order to perform 
an extrapolation N — > oo, we have analyzed the crossings 
between neighboring system sizes A2 > Ni > 9 x 9. This 
leads to the estimate 



T< 

C 

\J\ 



= 1.566 ±0.005 



(25) 



for the thermodynamic limit N — > 00. 



C. Nature of the phase transition for J < 

Having determined the transition temperature for J < 
0, one would like to clarify the universality class of the 
phase transition. 
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FIG. 14. (Color online) Probability to find a state with energy 
E/\J\ on the N = 90 x 90 lattice for J < at two selected 
temperatures: T/\J\ = 1.566 (left) and 1.7025 (right). Error 
bars in the lower panel do not exceed the width of the lines. 



On the one hand, there is no evidence for any latent 
heat in the specific heat at T<, see Fig. [lO^a), i.e., the 
ordering transition appears to be continuous for J < 0. 
On the other hand, a negative dip in the Binder cumu- 
lant for T > T c , as observed in Fig. [T3|[a) is sometimes 
taken as evidence for a first-order transition (see, e.g., 
Ref. EH) . In order to distinguish better between the two 
scenarios we use histograms of the energy E of the mi- 
crostates realized in the Monte-Carlo procedure P^HSZI We 
have collected such histograms for several system sizes 



and temperatures. Fig. 14 shows two representative cases 
on the N = 90 x 90 lattice, namely T = 1.7025 \J\ which 
corresponds to the maximum of the specific heat for the 
90x90 lattice (compare Fig.JlOja)) and T = 1.566 | J|, the 
estimated critical temperature of the infinite system, see 
Eq. (25). We always find bell-shaped almost Gaussian 
distributions, which are characteristic for a continuous 
transition. We never observed any signatures of a split- 
ting of this single peak into two, as would be expected 
for a first-order transition.ESHIU Hence, the transition ap- 
pears to be continuous and we will now try to charac- 
terize its universality class further in terms of critical 
exponents. 

We start with the correlation length exponent v which 
can be extracted from the finite-size behavior of the 
Binder cumulant: close to T c , the Binder cu mulant 
should scale with the linear size of the system L a d 46 ' 52 ' 53 1 
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FIG. 15. (Color online) Scaling of different quantities with 
linear size L for Lx L lattices, J < 0, and close to the critical 
temperature . The slope of the Binder cumulant yields 
the correlation length exponent v (top panel), the sublattice 
magnetization m s yields the exponent j3, and the specific heat 
C yields the exponent a. Lines show the fits which have been 
used to estimate the exponents. Note that the scale is double- 
logarithmic in all three panels. 



Fig. 13 b) shows that there is very little curvature in 
the Binder cumulants U s as a function of temperature T 
close to the estimated critical temperature ( 25 ) . There- 
fore dC/g/dT can be extracted without much sensitivity 
to the error of ( 25 ) . The result is shown by the symbols 

Now we can determine v by 



15 



in the top panel of Fig 
fitting this to ( 26 ) . Since inclusion of the correction term 
renders the fit 
(i.e., we set b - 



unstable, we use only the leading term 
in (26)). A fit (which is shown by the 
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line in the top panel of Fig. 15 1 then yields 



- = 0.24 ±0.02. 
v 



(27) 



We now turn to the order parameter exponent j3 which 
can be extracted from the finite-size behavior of the sub- 
lattice order parameter m s . The sublattice order param- 
eter should have a scaling behavior (see for example Refs. 
IMlandf 



(28) 



Specialization of ( 28 ) to T = T c yields 



T=T C 



L- 20 /»M 2 (O) 



(29) 



The middle panel of Fig. [15] shows the Monte-Carlo re- 
sults for m 2 at three temperatures which cover the es- 



timate (25) for Tf- and its error bars. The fits of these 



panel of Fig. 



results to ( 29 ) which are shown by the lines in the middle 



151 lead to 



2/3 



0.257 ± 0.006. 



(30) 



Finally, we turn to the specific heat exponent a. We 
proceed in the same manner as for the order parame- 
ter exponen t j3 an d make again a scaling ansatz for the 
specific heat 



C\ 



T=T a 



(31) 



The lower panel of Fig. 15 shows the specific heat results 
at the estimate ( 25 ) for T£* . One observes that this does 
not follow a power law very well. Indeed, it is known 
that non-scaling contributions to the specific heat can be 
important.^ However, including a constant in the ansatz 
(31) does not lead to a stable fit. We therefore fit only 
the data for L = 27, 36, 45, and 90 (lines in the lower 



panel of Fig. 15). This procedure leads to the estimate 



= 0.016 ±0.003. 



(32) 



Note that the error bar includes just the error of the fit. 
In view of the deviations from a simple power law, this 
is probably too optimistic. Indeed, a could very well 
be (slightly) negative, as is suggested by the fact that 
the maximal value of C in Fig. 10 'a) decreases when the 
system size increases from N = 45 x 45 to 90 x 90. 



Even if the error bars in (25), (27), and (32) should 



be too optimistic, it remains safe to conclude that we 
find a rather large correlation length exponent v > 3. 
It should be noted that in combination with a specific 
heat exponent a sa 0, we then find that the hyperscaling 



d v 



(33) 
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FIG. 16. (Color online) Average squared sublattice magne- 
tization for the classical model with J > 0. Error bars are 
of the order of the lines' width in this graph. Inset: aver- 
age squared sublattice magnetization in the low-temperature 
region. 



(with the spatial dimension d = 2) is strongly violated. 
On the other hand, we could use the relation (33) to 
estimate a, in particular if we expect it to be negative 
(compare, e.g., Ref. 1591 for a similar situation). Inser- 



tion of (27) into (33) yields a very negative exponent 
ajv = —1.52 ± 0.04. Again, this reflects the large ex- 
ponent v. In fact, a large correlation length exponent 
v 4 has be en fo und in other two-dimensional disor- 
dered systems! 60 * 61 ! However, in those cases the large ex- 
ponent corresponds to approaching the critical point via 
a fine-tuned direction in a two-dimensional parameter 
space and there is a se cond, substantially smaller cor- 
relation length exponent ! 60 * 61 ! Thus, we are left with not 
completely unreasonable, but definitely highly unusual 
values of the critical exponents v and a. 



D. Critical temperature for J > 

As mentioned earlier, for J > 0, the specific heat alone 
is only mildly conclusive regarding the existence of a 
low-temperature phase transition to a three-sublattice or- 
dered state. This statement requires to be supported by 
the analysis of other observables. The sublattice magne- 
tization ( 23 ) will tell us whether significant order is devel- 



oping in the low-temperature region or not. Our results 
shown in Fig. [16] show that three-sublattice order is in- 
deed developing, although an appreciable order develops 
only at temperatures that are so low that they become in- 
creasingly difficult to access with increasing system size. 
To take a closer look at the low-temperature ordered 
state, we took some snapshots of the system during the 
simulation for N — 12 x 12 spins. A typical configuration 
is reproduced in Fig. [8| While the global structure cor- 
responds indeed to a 120° three-sublattice ordered state, 
we also observe the presence of defects. The presence 
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FIG. 17. (Color online) Binder cumulant for the classical 
model with J > in the low-temperature region. Error bars 
are of the order of the lines' width in this graph. Inset: Binder 
cumulant for 6 x 6, 9 x 9, and 12 x 12 spins zoomed around 
the crossing region. 



of these defects is neither a surprise nor in contradiction 
with the existence of the phase transition, as they are in 
fact necessary ingredients of the order-by-disorder mech- 
anism. Note also that the spins in Fig. |8]lie essentially in 
the x — y-plane and -up to small fluctuations- are aligned 
with the lattice. 



As for J < 0, the Binder cumulant (24) allows us both 



to further support our conclusions concerning the low- 
temperature ordered state and to obtain an estimate of 
T>. Fig. 17 shows that all the curves for the different 
system sizes cross in a region around T rj 0.015J. First, 
this is a strong argument in favor of the existence of the 
ordering transition. We used the smallest three system 
sizes (N = 6 x 6, 9 x 9, and 12 x 12) for which we have 
the best statistics to obtain an estimate for the transition 
temperature: 



T> 

C 

J 



0.0125 ±0.0009. 



(34) 



The error bars on the data are unfortunately too large 
to get precise values for the critical exponents and thus 
prevent us from investigating the nature and the uni- 
versality class of the transition. However, the fact that 
the low-temperature ordered state breaks the same sym- 
metries irrespectively of the sign of J suggests that the 
universality class for J > is the same as for J < 0. 



VI. DISCUSSION AND CONCLUSIONS 

Although the Hamiltonian ([I]) may seem unusual in 
the context of frustrated magnetism, it is instructive in 
many respects and illustrates the rich phenomenology of- 
ten present in this subject. Either seen as the strong 
trimerization limit of the kagome lattice of spin 1/2 in 
a magnetic field, or as a possible illustration of an or- 



bital modelJ^HSZl the underlying physics associated to 
this Hamiltonian is extremely interesting for both signs 
of the coupling constant J. 

In the quantum case (for spin 1/2) we show, by study- 
ing the low-energy spectra using the Lanczos method, 
that a thermodynamic gap of the order of 3 | J| is present 
for J < 0, while for J > the gap, if present, would be 
at most of the order of 0.02 J. The six-fold degeneracy 
of the would-be ordered ground state which is predicted 
by semiclassical considerations is not observed in our nu- 
merical results, probably due to the small lattices con- 
sidered. These results illustrate very well how the deep 
quantum (S = 1/2) regime differs from the large S spin- 
wave predictions. The specific heat curves point to a 
phase transition around T w | J| for J < 0, while a lower 
temperature peak shows up in the positive J case. This 
last peak could be due to an ordering phase transition at 
a very low temperature T c < J/100. In both cases one 
is tempted to envisage a finite-temperature phase tran- 
sition whose nature could be understood by the analysis 
of the classical model. 

The analysis of the classical model has turned out to 
be also quite interesting and instructive. For J < 
the lowest-energy configuration consists in an in-plane 
antiferromagnetic arrangement of the spins with given 
chirality accompanied by a 'spurious' continuous rota- 
tional degeneracy which docs not correspond to any sym- 
metry of the Hamiltonian. This pseudo degeneracy is 
lifted by entropy at finite temperature giving rise to 
an ordering at low temperature as observed by Monte 
Carlo data which locate the transition temperature at 
T c /\ J\ = 1.566 ± 0.005. Inspection of the histograms of 
the energy close to the transition temperature gave no 
evidence of a first-order transition. Hence, we analyzed 
it within the scenario of a continuous transition and es- 
timated unusual values for the critical exponents a and 
v, strongly violating the hyperscaling relation. It should 
nevertheless be mentioned that we cannot exclude the ex- 
istence of a crossover scale which exceeds the lattice sizes 
accessible to us. The fact that lock-in of the spin com- 
ponents to the lattice requires a certain length scale may 
point in this direction. An unambiguous determination 
of the universality class of the transition would require 
improved methods. A first possibility is to restrict the 
degrees of freedom to the in-plane configurations^ which 
are realized in the low-temperature limit. Even more ef- 
ficiency could be gained by additionally restricting each 
spin variable to the 6 spin directions which are stabilized 
in the zero-temperature limit. However, it remains to 
be investigated whether the second modification changes 
the universality class of the transition. 

For J > the situation is even more interesting. Al- 
though a 'spurious' rotational degeneracy is also present 
for the antiferromagnetic 120° configuration (with the 
opposite chirality than the one for J < 0) , the manifold of 
lowest-energy configurations is more complex. There ex- 
ist local discrete 'flips' of triangles which bring one from 
the homogeneous antiferromagnetic lowest-energy config- 
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FIG. 18. (Color online) (a) Identification of chirality pseudo- 
spin states of a triangle with spin configurations on a triangle 
for the underlying kagome lattice, (b) Expression of the 120° 
ordered state on a triangle in the effective model in terms of 
spin configurations of the corresponding 9 sites in the under- 
lying kagome lattice. Note that chirality spins for the effec- 
tive model lie in the x-y chirality plane whereas spins on the 
kagome lattice point along the z-axis. 



tion of Fig. [TJb) for the chirality pseudo spins. Insertion 
into the 120° wavefunction for a triangle of the triangu- 
lar lattice on the left side of Fig. 18 'b) then yields the 
expression in terms of the 8 spin configurations of a nine- 
site unit of the underlying kagome lattice shown on the 
right side of Fig. 18 b). Note that the two terms on the 
first line of the right side of Fig. 18 'b) amount exactly to 



the variational wave function for the m agne tization 1/3 
state of the homogeneous kagome lattice} 7 -^ as it follows 
from a mapping to a quantum-dimer model on the honey- 
comb lattice™ Thus, the present results for the strongly 
trimerized kagome lattice may be smoothly connected to 
the 1/3 plateau state of the homogeneous kagome lattice. 

To conclude, the model ([I]) has turned out to be a 
very interesting laboratory to understand the emergence 
of a hierarchy of energy scales originating from differ- 
ent levels of order by disorder. The emergence of such 
a hierarchy in a classical model is related to similar hi- 
erarchies in quantum systems as for example the huge 
difference between the magnetic and non-magnetic gaps 
in the kagome spin 1/2 system at magnetization 1/3. 7 
Moreover, if the transitions observed in this work can be 
confirmed to be continuous, the exponents will probably 
correspond to exotic models, like parafermionic confor- 
mal field theories.^ 
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Appendix A: Relation to the Heisenberg model on 
the trimerized kagome lattice 



The Hamiltonian ( TH) has a lready been derived several 
times in the literature P 6 * 38 * 42 ! Nevertheless, for complete- 
ness we also give a derivation. 

We start from the interaction between the triangles of 
the spin-1/2 Heisenberg model on a trimerized kagome 
lattice 



^int — ^int <$A,i ' $B,j 



(Al) 



{id) 



— Jint I 9 (^A,i $B,j + $A,i ^B,j) + $A,i Sb.j \ i 

where the sum over runs over the nearest-neighbor 
pairs of triangles i and j in Fig. [lja) and the corners 
of the triangle A and B have to be chosen such as to 



17 



match the connecting bond. The 5a, i are physical spin- 
1/2 operators. 

In first order and for N t riangles, we need to compute 
the matrix elements of ( |Al[ ) between all 2 N combinations 
of the states in Fig. [ljb). Note that the expectation val- 
ues of the operators acting on different triangles factor- 
ize. Since in the present case magnetization is fixed in 
each triangle to 1/3, matrix elements of 5 A i vanish, thus 
simplifying the derivation considerably. 

For the lower left corner of a triangle i we find 



15, 



L.i 



\si ti \-) 



1/6 -oj 2 /3 

-w/3 1/6 

-(- 

3 V 2 



(A2) 



for the lower right corner 



I<SJU+) (+155,1- 



l-\sz 
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1/1_ T , 



3 12 



(A3) 



and finally for the top corner 



|5 

15? 



T,il 
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1 
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(A4) 



Using (A2)-(A4) for the matrix elements of ( Al ), we find 



the effective Hamiltonian 



H e ff. = 




-TV 



i-T' 
2 J 



E 

«fej» 



eg 



— Tif 



}_ _rpE 

2 3 



-2? 



, (A5) 



where the three sums run over the different bond direc- 
tions as sketched in Fig. [IT a). Since the sum over roots 
of unity vanishes, we have Tf + Tf + Tf = 0. Hence, 
Eq. ( A5 1 can be rewritten as 



H, 



eff. 




T A T B 



E 

[[Ml 



T C T B 



(A6) 



Up to an additive constant, this is nothing but Eq. ([I]) 
with J = Ji n t/9. The intra-triangle coupling needs to 
be chosen positive in order for the two states shown in 
Fig. |ljb) to be ground states of a triangle. Accordingly, 
it is natural to also choose the inter-triangle coupling J; nt 
positive, i.e., J > 0. 

Very similar arguments can be applied, e.g., to spin- 
less fermions with nearest-neighbor repulsionj^ leading 
to the same effective Hamiltonian. 
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